Vessel Densities in the Belgian Part of the North Sea
Lennert Schepers
4/28/2020
Vessel density in Belgium
In this small exercise we want to analyse vessel route densities in the Belgian part of the North Sea. The data for this exercise comes from EMODnet Human activities. EMODnet Human Activities has two types of vessel density data, one created by the Human Activities portal themselves, giving the vessel hours per square km per month by ship type. See here.
And one created by the European Maritime Safety Agency (EMSA), giving the number of routes per square km per month) by ship type. The advantage of the latter is that this provides recent information. For example, writing now 28th of April, the monthly aggregated data is already available for March 2020. For details, see here
First, we load the R libraries that we will use
## load rasters
library(raster)
library(sf)
library(mapview)
library(ggplot2)
library(data.table)For this exercise, we are only interested in Belgian waters, so we’ll use bounding box of the Belgian Exclusive Economic ZOne:
BCP <- st_read("http://geo.vliz.be/geoserver/MarineRegions/wfs?service=WFS&version=2.0.0&request=GetFeature&typeNames=eez&cql_filter=mrgid=%273293%27&outputFormat=application/json")## Reading layer `OGRGeoJSON' from data source `http://geo.vliz.be/geoserver/MarineRegions/wfs?service=WFS&version=2.0.0&request=GetFeature&typeNames=eez&cql_filter=mrgid=%273293%27&outputFormat=application/json' using driver `GeoJSON'
## Simple feature collection with 1 feature and 32 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2.238333 ymin: 51.08931 xmax: 3.370403 ymax: 51.87611
## geographic CRS: WGS 84# BCP <- st_cast(BCP, "POLYGON")
BCP_bbox <- st_bbox(BCP$geometry)
mapview(BCP)For the Scheldt Estuary:
WScheldt <- st_read("http://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.1.0&typename=MarineRegions:seavox_v16&outputformat=json&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Esub_region%3C%2FPropertyName%3E%3CLiteral%3EWESTERN+SCHELDT%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E")## Reading layer `OGRGeoJSON' from data source `http://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.1.0&typename=MarineRegions:seavox_v16&outputformat=json&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Esub_region%3C%2FPropertyName%3E%3CLiteral%3EWESTERN+SCHELDT%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E' using driver `GeoJSON'
## Simple feature collection with 1 feature and 23 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 3.546714 ymin: 51.32974 xmax: 4.252759 ymax: 51.46375
## geographic CRS: WGS 84WScheldt_bbox <- st_bbox(WScheldt$geometry)
mapview(WScheldt)bbox_BCP_Sch <- st_bbox(st_union(BCP,WScheldt))Now, we create a function to retrieve the route density for a particular fishing type, month and year in Belgian waters
base_url <- 'https://ows.emodnet-humanactivities.eu/wcs?SERVICE=WCS&VERSION=1.0.0&request=GetCapabilities'
# BACKGROUND INFO ON EMODNet Human Activities' WCS Service, and the 'route density maps'
# The Open Geospatial Consortium (OGC) Web Coverage Service Interface Standard (WCS) defines Web-based retrieval of coverages – that is, digital geospatial information representing space/time-varying phenomena.
#
# Get Capabilities:
# https://ows.emodnet-humanactivities.eu/wcs?SERVICE=WCS&VERSION=1.0.0&request=GetCapabilities
#
# Get Coverage Example:
# https://ows.emodnet-humanactivities.eu/wcs?service=wcs&version=1.0.0&request=getcoverage&coverage=emodnet:2019_01_rd_All&crs=EPSG:4326&BBOX=15,20.5,30,32.5&format=image/tiff&interpolation=nearest&resx=0.00833333&resy=0.00833333
#
# The names of the web service layers are -
# '<year>_<month>_rd_<vessel code>' Eg. '2019_01_rd_All' for all vessels during January 2019
# '<year>_<season>_rd_<vessel code>' Eg. '2019_spring_rd_All' for all vessels during Spring 2019
# '<year>_yearly_rd_<vessel code>' Eg. '2019_yearly_rd_All' for the total of all vessels during 2019.
# The vessel codes are as follows:
#
# All - All types
# 01 - Cargo
# 02 - Fishing
# 03 - Passenger
# 04 - Tanker
# 05 - Other
#
#Available data: January 2019 - March 2020
# function to get belgian raster data, based on year, month, vessel_code
getBelgianCSW <- function(year, month, vessel_code){
url <- paste0("https://ows.emodnet-humanactivities.eu/wcs?service=wcs&version=1.0.0&request=getcoverage&coverage=emodnet:",
year, "_",
formatC(month, width = 2, flag = 0), "_rd_",
vessel_code, "&",
"crs=EPSG:4326&BBOX=", paste(bbox_BCP_Sch, collapse = ","), "&",
"format=image/tiff&",
"interpolation=nearest&resx=0.00833333&resy=0.00833333"
)
raster(url)
}Now we extract the Vessel density data from the WCS function, and store it in a rasterbrick:
Cargo <- brick(x = c('Cargo-2019-01' = getBelgianCSW(2019, 1, "01"),
'Cargo-2019-02' = getBelgianCSW(2019, 2, "01"),
'Cargo-2019-03' = getBelgianCSW(2019, 3, "01"),
'Cargo-2019-04' = getBelgianCSW(2019, 4, "01"),
'Cargo-2019-05' = getBelgianCSW(2019, 5, "01"),
'Cargo-2019-06' = getBelgianCSW(2019, 6, "01"),
'Cargo-2019-07' = getBelgianCSW(2019, 7, "01"),
'Cargo-2019-08' = getBelgianCSW(2019, 8, "01"),
'Cargo-2019-09' = getBelgianCSW(2019, 9, "01"),
'Cargo-2019-10' = getBelgianCSW(2019, 10, "01"),
'Cargo-2019-11' = getBelgianCSW(2019, 11, "01"),
'Cargo-2019-12' = getBelgianCSW(2019, 12, "01"),
'Cargo-2020-01' = getBelgianCSW(2020, 1, "01"),
'Cargo-2020-02' = getBelgianCSW(2020, 2, "01"),
'Cargo-2020-03' = getBelgianCSW(2020, 3, "01"),
'Cargo-2020-04' = getBelgianCSW(2020, 4, "01")
)
)
Fishing <- brick(x = c('Fishing-2019-01' = getBelgianCSW(2019, 1, "02"),
'Fishing-2019-02' = getBelgianCSW(2019, 2, "02"),
'Fishing-2019-03' = getBelgianCSW(2019, 3, "02"),
'Fishing-2019-04' = getBelgianCSW(2019, 4, "02"),
'Fishing-2019-05' = getBelgianCSW(2019, 5, "02"),
'Fishing-2019-06' = getBelgianCSW(2019, 6, "02"),
'Fishing-2019-07' = getBelgianCSW(2019, 7, "02"),
'Fishing-2019-08' = getBelgianCSW(2019, 8, "02"),
'Fishing-2019-09' = getBelgianCSW(2019, 9, "02"),
'Fishing-2019-10' = getBelgianCSW(2019, 10, "02"),
'Fishing-2019-11' = getBelgianCSW(2019, 11, "02"),
'Fishing-2019-12' = getBelgianCSW(2019, 12, "02"),
'Fishing-2020-01' = getBelgianCSW(2020, 1, "02"),
'Fishing-2020-02' = getBelgianCSW(2020, 2, "02"),
'Fishing-2020-03' = getBelgianCSW(2020, 3, "02"),
'Fishing-2020-04' = getBelgianCSW(2020, 4, "02")
)
)
mapview(Fishing[[1]])Passenger <- brick(x = c('Passenger-2019-01' = getBelgianCSW(2019, 1, "03"),
'Passenger-2019-02' = getBelgianCSW(2019, 2, "03"),
'Passenger-2019-03' = getBelgianCSW(2019, 3, "03"),
'Passenger-2019-04' = getBelgianCSW(2019, 4, "03"),
'Passenger-2019-05' = getBelgianCSW(2019, 5, "03"),
'Passenger-2019-06' = getBelgianCSW(2019, 6, "03"),
'Passenger-2019-07' = getBelgianCSW(2019, 7, "03"),
'Passenger-2019-08' = getBelgianCSW(2019, 8, "03"),
'Passenger-2019-09' = getBelgianCSW(2019, 9, "03"),
'Passenger-2019-10' = getBelgianCSW(2019, 10, "03"),
'Passenger-2019-11' = getBelgianCSW(2019, 11, "03"),
'Passenger-2019-12' = getBelgianCSW(2019, 12, "03"),
'Passenger-2020-01' = getBelgianCSW(2020, 1, "03"),
'Passenger-2020-02' = getBelgianCSW(2020, 2, "03"),
'Passenger-2020-03' = getBelgianCSW(2020, 3, "03"),
'Passenger-2020-04' = getBelgianCSW(2020, 4, "03")
)
)
Tanker <- brick(x = c('Tanker-2019-01' = getBelgianCSW(2019, 1, "04"),
'Tanker-2019-02' = getBelgianCSW(2019, 2, "04"),
'Tanker-2019-03' = getBelgianCSW(2019, 3, "04"),
'Tanker-2019-04' = getBelgianCSW(2019, 4, "04"),
'Tanker-2019-05' = getBelgianCSW(2019, 5, "04"),
'Tanker-2019-06' = getBelgianCSW(2019, 6, "04"),
'Tanker-2019-07' = getBelgianCSW(2019, 7, "04"),
'Tanker-2019-08' = getBelgianCSW(2019, 8, "04"),
'Tanker-2019-09' = getBelgianCSW(2019, 9, "04"),
'Tanker-2019-10' = getBelgianCSW(2019, 10, "04"),
'Tanker-2019-11' = getBelgianCSW(2019, 11, "04"),
'Tanker-2019-12' = getBelgianCSW(2019, 12, "04"),
'Tanker-2020-01' = getBelgianCSW(2020, 1, "04"),
'Tanker-2020-02' = getBelgianCSW(2020, 2, "04"),
'Tanker-2020-03' = getBelgianCSW(2020, 3, "04"),
'Tanker-2020-04' = getBelgianCSW(2020, 4, "04")
)
)
Other <- brick(x = c('Other-2019-01' = getBelgianCSW(2019, 1, "05"),
'Other-2019-02' = getBelgianCSW(2019, 2, "05"),
'Other-2019-03' = getBelgianCSW(2019, 3, "05"),
'Other-2019-04' = getBelgianCSW(2019, 4, "05"),
'Other-2019-05' = getBelgianCSW(2019, 5, "05"),
'Other-2019-06' = getBelgianCSW(2019, 6, "05"),
'Other-2019-07' = getBelgianCSW(2019, 7, "05"),
'Other-2019-08' = getBelgianCSW(2019, 8, "05"),
'Other-2019-09' = getBelgianCSW(2019, 9, "05"),
'Other-2019-10' = getBelgianCSW(2019, 10, "05"),
'Other-2019-11' = getBelgianCSW(2019, 11, "05"),
'Other-2019-12' = getBelgianCSW(2019, 12, "05"),
'Other-2020-01' = getBelgianCSW(2020, 1, "05"),
'Other-2020-02' = getBelgianCSW(2020, 2, "05"),
'Other-2020-03' = getBelgianCSW(2020, 3, "05"),
'Other-2020-04' = getBelgianCSW(2020, 4, "05")
)
)
All <- brick(x = c('All-2019-01' = getBelgianCSW(2019, 1, "All"),
'All-2019-02' = getBelgianCSW(2019, 2, "All"),
'All-2019-03' = getBelgianCSW(2019, 3, "All"),
'All-2019-04' = getBelgianCSW(2019, 4, "All"),
'All-2019-05' = getBelgianCSW(2019, 5, "All"),
'All-2019-06' = getBelgianCSW(2019, 6, "All"),
'All-2019-07' = getBelgianCSW(2019, 7, "All"),
'All-2019-08' = getBelgianCSW(2019, 8, "All"),
'All-2019-09' = getBelgianCSW(2019, 9, "All"),
'All-2019-10' = getBelgianCSW(2019, 10, "All"),
'All-2019-11' = getBelgianCSW(2019, 11, "All"),
'All-2019-12' = getBelgianCSW(2019, 12, "All"),
'All-2020-01' = getBelgianCSW(2020, 1, "All"),
'All-2020-02' = getBelgianCSW(2020, 2, "All"),
'All-2020-03' = getBelgianCSW(2020, 3, "All"),
'All-2020-04' = getBelgianCSW(2020, 4, "All")
)
)No we create a data frame to store the statistics that we will extract from the rasters. We extract for all values of the Belgian EEZ, the average value of all raster cells.
# create data frame, with 1st column 'time'
df.sr <- data.frame(Date = seq(as.Date("2019/1/1"), by = "month", length.out = 16))
# extract all values from the rasters, by the BCP polygon, and average the values.
df.sr$Cargo <- as.vector(extract(Cargo, BCP, fun = mean))
df.sr$Fishing <- as.vector(extract(Fishing, BCP, fun = mean))
df.sr$Passenger <- as.vector(extract(Passenger, BCP, fun = mean))
df.sr$Tanker <- as.vector(extract(Tanker, BCP, fun = mean))
df.sr$Other <- as.vector(extract(Other, BCP, fun = mean))
df.sr$All <- as.vector(extract(All, BCP, fun = mean))
# transform to long format:
df.sr.long <- melt.data.table(as.data.table(df.sr),
id.vars = "Date")A plot of the result:
ggplot(df.sr.long, aes(x = Date, y = value, group = variable)) +
geom_point(aes(color = variable)) +
geom_line(aes(color = variable)) +
labs(title = "Vessel route densities in the Belgian part of the North Sea",
subtitle = "Monthly average route density (routes/km²/ship type) for the Belgian EEZ",
x = "Month",
y = "Average routes / km² / ship type") +
theme_minimal()
Scheldt estuary:
# create data frame, with 1st column 'time'
df.Scheldt <- data.frame(Date = seq(as.Date("2019/1/1"), by = "month", length.out = 16))
# extract all values from the rasters, by the WScheldt polygon, and average the values.
df.Scheldt$Cargo <- as.vector(extract(Cargo, WScheldt, fun = mean))
df.Scheldt$Fishing <- as.vector(extract(Fishing, WScheldt, fun = mean))
df.Scheldt$Passenger <- as.vector(extract(Passenger, WScheldt, fun = mean))
df.Scheldt$Tanker <- as.vector(extract(Tanker, WScheldt, fun = mean))
df.Scheldt$Other <- as.vector(extract(Other, WScheldt, fun = mean))
df.Scheldt$All <- as.vector(extract(All, WScheldt, fun = mean))
# transform to long format:
df.Scheldt.long <- melt.data.table(as.data.table(df.Scheldt),
id.vars = "Date")A plot of the result:
ggplot(df.Scheldt.long, aes(x = Date, y = value, group = variable)) +
geom_point(aes(color = variable)) +
geom_line(aes(color = variable)) +
labs(title = "Vessel route densities in the Western Scheldt",
subtitle = "Monthly average route density (routes/km²/ship type) in the Western Scheldt",
x = "Month",
y = "Average routes / km² / ship type") +
theme_minimal()
Density Maps
If we want to compare the difference between 2020 vs 2019:
# Average map of February-April 2019
r2019_0204 <- brick(c(
Cargo_2019_0204 = mean(Cargo[[2:4]]),
Fishing_2019_0204 = mean(Fishing[[2:4]]),
Passenger_2019_0204 = mean(Passenger[[2:4]]),
Tanker_2019_0204 = mean(Tanker[[2:4]]),
Other_2019_0204 = mean(Other[[2:4]]),
All_2019_0204 = mean(All[[2:4]])
)
)
# Average map of February-April 2020
r2020_0204 <- brick(c(
Cargo_2020_0204 = mean(Cargo[[14:16]]),
Fishing_2020_0204 = mean(Fishing[[14:16]]),
Passenger_2020_0204 = mean(Passenger[[14:16]]),
Tanker_2020_0204 = mean(Tanker[[14:16]]),
Other_2020_0204 = mean(Other[[14:16]]),
All_2020_0204 = mean(All[[14:16]])
)
)
mapview(r2019_0204) + mapview(r2020_0204)Anomaly maps
If we want to compare the difference between 2020 vs 2019:
# Anomaly map January-April 2020 - January-April 2019
Anomaly_Cargo <- brick(c(
Anomaly_Cargo_01 = mean(Cargo[[1]] - Cargo[[13]]),
Anomaly_Cargo_02 = mean(Cargo[[2]] - Cargo[[14]]),
Anomaly_Cargo_03 = mean(Cargo[[3]] - Cargo[[15]]),
Anomaly_Cargo_04 = mean(Cargo[[4]] - Cargo[[16]])
)
)
Anomaly_Fishing <- brick(c(
Anomaly_Fishing_01 = mean(Fishing[[1]] - Fishing[[13]]),
Anomaly_Fishing_02 = mean(Fishing[[2]] - Fishing[[14]]),
Anomaly_Fishing_03 = mean(Fishing[[3]] - Fishing[[15]]),
Anomaly_Fishing_03 = mean(Fishing[[4]] - Fishing[[16]])
)
)
Anomaly_Passenger <- brick(c(
Anomaly_Passenger_01 = mean(Passenger[[1]] - Passenger[[13]]),
Anomaly_Passenger_02 = mean(Passenger[[2]] - Passenger[[14]]),
Anomaly_Passenger_03 = mean(Passenger[[3]] - Passenger[[15]]),
Anomaly_Passenger_04 = mean(Passenger[[4]] - Passenger[[16]])
)
)
Anomaly_Tanker <- brick(c(
Anomaly_Tanker_01 = mean(Tanker[[1]] - Tanker[[13]]),
Anomaly_Tanker_02 = mean(Tanker[[2]] - Tanker[[14]]),
Anomaly_Tanker_03 = mean(Tanker[[3]] - Tanker[[15]]),
Anomaly_Tanker_04 = mean(Tanker[[4]] - Tanker[[16]])
)
)
Anomaly_Other <- brick(c(
Anomaly_Other_01 = mean(Other[[1]] - Other[[13]]),
Anomaly_Other_02 = mean(Other[[2]] - Other[[14]]),
Anomaly_Other_03 = mean(Other[[3]] - Other[[15]]),
Anomaly_Other_04 = mean(Other[[4]] - Other[[16]])
)
)
Anomaly_All <- brick(c(
Anomaly_All_01 = mean(All[[1]] - All[[13]]),
Anomaly_All_02 = mean(All[[2]] - All[[14]]),
Anomaly_All_03 = mean(All[[3]] - All[[15]]),
Anomaly_All_04 = mean(All[[4]] - All[[16]])
)
)
mapview(Anomaly_All)Plot anomalies January-March combined:
mapview(r2020_0204 - r2019_0204)# Anomaly graph?
#### TO DO #####